Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)Study 1
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))
df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))load("../data/Gaggero2021/data.rda")
df2 <- data |>
mutate(Gender = as.character(gender)) |>
mutate(across(starts_with("IAS "), as.numeric)) |>
rename(
Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
rename(
Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
)df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
rename(
Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
)df5_attention <- readxl::read_excel("../data/Arslanova2022/Participants2.xlsx",
sheet = "FINAL") |>
filter(ActorHR == 1) |>
select(Subj.ID, Gender, Age) |>
mutate(Gender = as.numeric(Gender), Age = as.numeric(Age)) |>
filter(!is.na(Age)) |> #error on the documentation of this participant - attention check failed but not noted
mutate(Gender = case_when(
Gender== 0 ~ "Male",
Gender== 1 ~ "Female" # based on paper reporting 65 women
))New names:
• `Why No?` -> `Why No?...4`
• `Why No?` -> `Why No?...6`
• `` -> `...12`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`
df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx") |>
filter(str_detect(Question.Key, "IAC_")) |>
filter(str_detect(Question.Key, "-quantised")) |>
pivot_wider(names_from = Question.Key, values_from = Response) |>
rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |>
rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |>
rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast,
Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar,
Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |>
rename(Subj.ID = "Participant.Public.ID") |>
select(-starts_with("Participant"))
df5 <- df5 |>
filter(Subj.ID %in% df5_attention$Subj.ID) |>
select(-Subj.ID) |>
mutate(across(everything(), as.numeric))df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |>
select(-SERIAL) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |>
select(-COUNTRY_OTHER) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)
df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |>
select(-COUNTRY_OTHER, -EDUCATION_OTHER, -Sample) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |>
select(-ID) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
)df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
select(-sex) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
rename(
Age = age,
Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
Affective_touch = Touch, Itch = ITCH
)
df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
rename(
Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
Affective_touch = Touch, Itch = ITCH)df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
rename(
Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
) |>
filter(!is.na(Urinate))
df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
rename(
Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
typeSample = Sample
)
df10c <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionPrimals/refs/heads/main/data/data_participants.csv") |>
select("Participant" = "participant_id", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI"))) |>
rename(
Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
) |>
filter(!is.na(Heart)) # participant is outlier and had total IAS scores below 0.4
df10d <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/rawdata_participants.csv") |>
select("Participant", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI", "PHQ4"))) |>
rename(
Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
) |>
filter(Participant %in% read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/data_participants.csv")$Participant)# pi_vars <- c("PI_Enticing", "PI_Alive", "PI_Safe", "PI_Good", "PI_Changing",
# "PI_Hierarchical", "PI_Understandable", "PI_AboutMe", "PI_Abundant",
# "PI_Acceptable", "PI_Beautiful", "PI_Cooperative", "PI_Funny",
# "PI_Harmless", "PI_Improvable", "PI_Intentional", "PI_Interconnected",
# "PI_Interesting", "PI_Just", "PI_Meaningful", "PI_NeedsMe",
# "PI_Pleasurable", "PI_Progressing", "PI_Regenerative", "PI_Stable",
# "PI_WorthExploring")
#
#
# correlation::correlation(
# df10c[, "IAS_Total", drop = FALSE],
# data2 = select(df10c, all_of(pi_vars)), p_adjust = "none"
# ) |>
# # Formatting correlation results
# mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
#
# # Releveling factors for visualization
# mutate(Parameter2 = fct_relevel(Parameter2, rev(pi_vars)),
# Parameter1 = fct_relevel(Parameter1, "IAS_Total")) |>
#
# # Plotting the correlation matrix
# ggplot(aes(x = Parameter1, y = Parameter2)) +
# geom_tile(aes(fill = r), color = "white") +
# geom_text(aes(label = val), size = 3) +
# labs(title = "Correlation Matrix") +
# scale_fill_gradient2(low = "blue", high = "red", mid = "white", limit = c(-1, 1), midpoint = 0, guide = guide_colourbar(ticks = FALSE)) +
# theme_minimal() +
# theme(
# legend.title = element_blank(),
# axis.title.x = element_blank(),
# axis.title.y = element_blank(),
# axis.text.x = element_text(hjust = 1)
# )df11a <- haven::read_sav("../data/Giulia/Interoceptive Attention ESM.sav") |>
select(AGE, GENDER, contains("IAS_ACC"), -IAS_ACC_20, starts_with(c("MAIA", "TAS20", "DEP", "ANX")), -ANX_21, -TAS20_21) |> #anx_21 , tas20_21 are attention checks
mutate_all(as.numeric) |>
rename(
Age = AGE, Gender = GENDER,
Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4,
Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7, Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9,
Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18,
Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_21, Itch = IAS_ACC_22
) |>
mutate(Gender = case_when(
Gender == 1 ~ "Male",
Gender == 2 ~ "Female",
Gender %in% c(3, 5) ~ "Other",
Gender == 4 ~ "NA"
)) |> ## compute dimensions for the TAS
mutate(TAS_DIF = TAS20_1 + TAS20_3 + TAS20_6 + TAS20_7 + TAS20_9 + TAS20_13 + TAS20_14,
TAS_DDF = TAS20_2 + (1-TAS20_4) + TAS20_11 + TAS20_12 + TAS20_17,
TAS_EOT = (1-TAS20_5) + TAS20_8 + (1-TAS20_10) + TAS20_15 + TAS20_16 + (1-TAS20_18) + (1-TAS20_19) + TAS20_20)
df11b <- haven::read_sav("../data/Giulia/Sleep and Exteroceptive Interoceptive Sensitivity.sav") |> # 165 participants
filter(Total_Attention_Fail == 0) |> # 32 participants removed based on failing one or more checks
select(AGE, GENDER, contains("IAS_ACC"), ) |>
mutate_all(as.numeric) |>
rename( Age = AGE, Gender = GENDER,
Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4, Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7,
Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9, Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18, Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_20,
Itch = IAS_ACC_21 ) |>
mutate(Gender = case_when(
Gender == 1 ~ "Male",
Gender == 2 ~"Female",
Gender %in% c(3, 5) ~ "Other",
Gender == 4 ~ "NA", )) |>
na.omit()
# Save all
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b, df10c=df10c, df11a=df11a, df11b=df11b)
save(data, file = "../data/data.RData")Total N = 32256.
desc <- function(df){
summary <- df |> summarise(
n = n(),
mean = mean(Age, na.rm = TRUE),
sd = sd(Age, na.rm = TRUE),
min = min(Age, na.rm = TRUE),
max = max(Age, na.rm = TRUE),
female = sum(Gender == "Female", na.rm = TRUE)
)
return(summary)
}
desc_df1a <- desc(df1a)
desc_df1b <- desc(df1b)
desc_df2 <- desc(df2)
desc_df3 <- desc(df3)
#desc_df4 <- desc(df4) # no demographic data available
desc_df5 <- desc(df5_attention)
desc_df6 <- desc(df6)
desc_df7a <- desc(df7a)
desc_df7b <- desc(df7b)
desc_df7c <- desc(df7c)
desc_df8a <- desc(df8a)
desc_df8b <- desc(df8b)
desc_df9 <- desc(df9)
desc_df10a <- desc(df10a)
desc_df10b <- desc(df10b)
desc_df10c <- desc(df10c)
desc_df11a <- desc(df11a)
desc_df11b <- desc(df11b)
desc_total <- rbind(desc_df1a, desc_df1b, desc_df2, desc_df3, desc_df5, desc_df6, desc_df7a, desc_df7b, desc_df7c, desc_df8a, desc_df8b, desc_df9, desc_df10a, desc_df10b, desc_df10c, desc_df11a, desc_df11b)
## calculate weighted mean
total_mean = sum(desc_total$mean * desc_total$n) / sum(desc_total$n)
total_sd = sum(desc_total$sd * desc_total$n) / sum(desc_total$n)
perce_female = sum(desc_total$female)/sum(desc_total$n) *100library(gt)
# APA style ####
gt_apastyle <- function(gt_table, font.size=12) {
gt_table |>
gt::opt_table_lines(extent = "none") %>%
gt::tab_options(
heading.border.bottom.width = 2,
heading.border.bottom.color = "black",
heading.border.bottom.style = "solid",
table.border.top.color = "black",
table.border.top.style = "solid",
table.border.top.width = 2,
table_body.hlines.color = "white",
table_body.border.top.color = "black",
table_body.border.top.style = "solid",
table_body.border.top.width = 2,
heading.title.font.size = font.size,
table.font.size = font.size,
heading.subtitle.font.size = font.size,
table_body.border.bottom.color = "black",
table_body.border.bottom.width = 2,
table_body.border.bottom.style = "solid",
column_labels.border.bottom.color = "black",
column_labels.border.bottom.style = "solid",
column_labels.border.bottom.width = 1,
latex.use_longtable = FALSE
) |>
gt::opt_table_font(font = "times")
}
make_age <- function(age) {
age <- as.numeric(age)
mean(age, na.rm = TRUE) |>
insight::format_value(digits=1) |>
paste0(" ± ",
insight::format_value(sd(age, na.rm = TRUE), digits=1))
}
# Create the dataframe
table <- data.frame(
Reference = c('Murphy et al., (2020)', '', '', 'Gaggero et al., (2021)', 'Campos et al., (2022)', 'Todd et al., (2022)', 'Arslanova et al., (2022)', 'Brand et al., (2022)', 'Brand et al., (2023)', '', '', '', 'Lin et al., (2023)', '', '', 'VonMohr et al., (2023)', 'Makowski et al., (2023a)', 'Makowski et al., (2023b)', 'Makowski et al., (2023c)', 'Poerio et al., (2024)', 'Poerio et al., unpublished', 'Total'),
Sample = c('', 'Sample 1a', 'Sample 1b', 'Sample 2', 'Sample 3', 'Sample 4', 'Sample 5', 'Sample 6', '', 'Sample 7a', 'Sample 7b', 'Sample 7c', '', 'Sample 8a', 'Sample 8b', 'Sample 9', 'Sample 10', 'Sample 11', 'Sample 12', 'Sample 13', 'Sample 14', 'Sample 15'),
Language = c('', 'English', 'English' , 'English and Italian', 'Portuguese', 'English', 'English', 'German', '', 'German', 'German', 'German', '', 'Chinese', 'Chinese', 'English', 'English', 'English', 'English', 'English', 'English',''),
N = c('', nrow(df1a), nrow(df1b), nrow(df2), nrow(df3), nrow(df4), nrow(df5_attention), nrow(df6),'', nrow(df7a), nrow(df7b), 802,'', nrow(df8a), nrow(df8b), nrow(df9), nrow(df10a), nrow(df10b), nrow(df10c), nrow(df11a), nrow(df11b), 32214),
'Difference' = c('','', '', '', '', '', '', '', '', '', '', '', '', 'Collapsed "Itch" and "Tickling"', 'Collapsed "Itch" and "Tickling"', '', 'Analog scales', 'Analog scales', 'Analog scales', '', '', ''),
Age = c('', make_age(df1a$Age), make_age(df1b$Age), make_age(df2$Age), make_age(df3$Age), "48.6.6 ± 14.1*", make_age(df5_attention$Age), make_age(df6$Age),'', make_age(df7a$Age), make_age(df7b$Age), make_age(df7c$Age),'', make_age(df8a$Age), make_age(df8b$Age), make_age(df9$Age), make_age(df10a$Age), make_age(df10b$Age), make_age(df10c$Age), make_age(df11a$Age), make_age(df11b$Age), '48.6 ± 13.1'),
Range = c('', '18-69', '18-91', '18-58', '18-72', '18-92*', '18-73', '18-78', '', '18-79', '16-81', '18-72','', '16-60', '20-60', '18-93', '18-73', '17-76','18-50', '18-57', '18-60', '17-93'),
Female_Percentage = c('', '69.4%', '70.1%', '60.3%', '59.6%', '50%*', '46.8%', '78.7%', '', '79.5%', '77.7%', '68.9%', '', '57.0%', '56.2%', '73.2%', '50.3%', '53.0%', '76%', '74.8%', '75.9%', '71.6%'),
Availability= c('osf.io/3m5nh', '', '', 'osf.io/5x9sg', 'osf.io/j6ef3', 'osf.io/ms354', 'osf.io/mp3cy', 'osf.io/xwz6g', 'osf.io/3f2h6', '', '', '','osf.io/3eztd','', '', 'osf.io/7p9u5', 'github.com/RealityBending/IllusionGameReliability', 'github.com/DominiqueMakowski/PHQ4R', 'github.com/RealityBending/InteroceptionPrimals', 'osf.io/49wbv', '','')
)
table_apa <- table |>
gt() |>
cols_align(align = c("right"), columns = "Age") |>
cols_label(Age = "Age (Mean ± SD)", Female_Percentage = "Female %") |>
tab_footnote("* Information taken from the sample description of relevant paper rather than recomputed.")
gt_apastyle(table_apa, font.size=12)
gtsave(gt_apastyle(table_apa, font.size=9), "figures/table1.tex")
# gtsave(table_apa, "figures/table1.png")The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.
vars <- c(
"Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
"Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)
dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
mutate(Sample = "Sample 10a",
x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 10b",
x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10c <- estimate_density(select(df10c, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 10c",
x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens11a <- estimate_density(select(df11a, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 11a")
dens11b <- estimate_density(select(df11b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 11b")
p1_data <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b, dens10c, dens11a, dens11b)
linetype <- setNames(rep("solid", length(vars)), vars)
linetype["Affective touch"] <- "dotted"
linetype["Blood sugar"] <- "dashed"
linetype["Sex arousal"] <- "solid"
linetype["Bruise"] <- "dashed"
# proportion of non-1 and non-5 values
p1_data |>
filter(x > 1 & x < 5) |> # Keep only x values strictly between 1 and 5
summarise(proportion = sum(y) / sum(p1_data$y)) # Normalize by total density proportion
1 0.9994812
p1 <- p1_data |>
mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b", "Sample 10c", "Sample 11a", "Sample 11b"),
Parameter = ifelse(Parameter == "Affective_touch", "Affective touch", Parameter),
Parameter = ifelse(Parameter == "Blood_Sugar", "Blood sugar", Parameter),
Parameter = ifelse(Parameter == "Sex_arousal", "Sex arousal", Parameter)) |>
# mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |>
ggplot(aes(x = x, y = y, color = Parameter)) +
geom_line(aes(linetype = Parameter), linewidth = 0.7) +
scale_color_material_d() +
scale_linetype_manual(values = linetype) +
labs(x = "Response", y = "Distribution", color = "Item", linetype = "Item", title = "Item Distribution") +
guides(color = guide_legend(ncol = 1)) +
facet_wrap(~Sample, scales = "free_y", nrow = 3) +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(face="bold"))
p1data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))
data10c <- select(df10c, all_of(dens10c$Parameter))
data11a <- normalize(select(df11a, all_of(dens11a$Parameter)), range = c(1, 5))
data11b <- na.omit(normalize(select(df11b, all_of(dens11b$Parameter)), range = c(1, 5))) # remove NAS
row.names(data11b) <- NULL
data_all <- rbind(
data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
mutate(data8a, Tickle = NA),
mutate(data8b, Tickle = NA), data9,
mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA),
data10b, data10c, data11a, data11b
) An overall positive intercorrelation pattern, with no clear structure emerging.
make_correlation <- function(df) {
correlation::correlation(df, redundant = TRUE) |>
correlation::cor_sort() |>
# correlation::cor_lower() |>
mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE)),
r = ifelse(Parameter1 == Parameter2, NA, r),
lbl = ifelse(Parameter1 == Parameter2, "", format_value(r, zap_small = TRUE, lead_zero = FALSE)),
lbl = ifelse(p > .05, "", lbl),
Param1 = fct_relevel(str_replace(Parameter1, "_", " "), str_replace(levels(Parameter1), "_", " ")),
Param2 = fct_relevel(str_replace(Parameter2, "_", " "), str_replace(levels(Parameter2), "_", " "))) |>
# mutate(Parameter2 = fct_rev(Parameter2)) |>
ggplot(aes(x = Param1, y = Param2)) +
geom_tile(aes(fill = r), color = "white") +
geom_text(aes(label = lbl), size = 3) +
labs(title = "Correlation Matrix", subtitle = paste0("N = ", nrow(df))) +
scale_fill_gradientn(
colors = c("white", "#FFF9C4", "#FFF59D", "#FFEB3B", "#FFCA28", "#FF9800", "#FF5722", "#F44336", "#E91E63", "#9C27B0", "#673AB7", "#512DA8", "#311B92"),
na.value = "#311B92",
limits = c(0, 0.8),
guide = guide_colorbar(ticks.colour = NA)
) +
# scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
theme_minimal() +
theme(
legend.title = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold")
)
}
p2 <- make_correlation(data_all)
p2fig1 <- p1 / p2
ggsave("figures/Figure1.png", fig1, width=10, height=12, dpi=300, bg="white")make_correlation(data1a)make_correlation(data1b)make_correlation(data2)make_correlation(data3)make_correlation(data4)make_correlation(data5)make_correlation(data6)make_correlation(data7a)make_correlation(data7b)make_correlation(data7c)make_correlation(data8a)make_correlation(data8b)make_correlation(data9)make_correlation(data10a)make_correlation(data10b)make_correlation(data10c)make_correlation(data11a)make_correlation(data11b)Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).
uva0 <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva0Variable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.363
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.290
Urinate Defecate 0.265
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Sneeze Cough 0.230
Heart Breathing 0.225
Hungry Thirsty 0.217
uva0$keep_remove$keep
[1] "Itch"
$remove
[1] "Tickle"
uva1a <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva1aVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Tickle Itch 0.286
Wind Burp 0.275
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Sneeze Cough 0.244
Urinate Defecate 0.241
uva1a$keep_removeNULL
uva1b <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva1bVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.350
Sneeze Cough 0.317
Wind Burp 0.309
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Urinate Defecate 0.278
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva1b$keep_remove$keep
[1] "Cough" "Burp" "Itch"
$remove
[1] "Sneeze" "Wind" "Tickle"
uva2 <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva2Variable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Tickle Itch 0.283
Heart Breathing 0.252
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Sneeze Cough 0.230
Urinate Defecate 0.223
Hungry Thirsty 0.200
uva2$keep_removeNULL
uva3 <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva3Variable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.445
Sneeze Cough 0.318
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.256
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Bruise Blood_Sugar 0.219
Urinate Defecate 0.217
Heart Breathing 0.208
uva3$keep_remove$keep
[1] "Sneeze" "Itch"
$remove
[1] "Cough" "Tickle"
uva4 <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva4Variable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Heart Breathing 0.415
Tickle Itch 0.351
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Hungry Thirsty 0.293
Urinate Defecate 0.282
Wind Burp 0.275
Sneeze Cough 0.251
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Muscles Pain 0.205
uva4$keep_remove$keep
[1] "Heart" "Itch"
$remove
[1] "Breathing" "Tickle"
uva5 <- EGAnet::UVA(data = data5, cut.off = 0.3)Warning: Some variables did not belong to a dimension: Blood_Sugar, Affective_touch
Use caution: These variables have been removed from the TEFI calculation
uva5Variable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Tickle Itch 0.247
Bruise Itch 0.246
Breathing Vomit 0.207
uva5$keep_removeNULL
uva6 <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva6Variable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Urinate Defecate 0.416
Sneeze Cough 0.332
Wind Burp 0.314
Tickle Itch 0.303
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Hungry Thirsty 0.274
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Heart Breathing 0.237
uva6$keep_remove$keep
[1] "Urinate" "Sneeze" "Wind" "Itch"
$remove
[1] "Defecate" "Cough" "Burp" "Tickle"
uva7a <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva7aVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Sneeze Cough 0.434
Tickle Itch 0.376
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.295
Urinate Defecate 0.294
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva7a$keep_remove$keep
[1] "Sneeze" "Tickle"
$remove
[1] "Cough" "Itch"
uva7b <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva7bVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Urinate Defecate 0.317
Tickle Itch 0.311
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Heart Breathing 0.281
Sneeze Cough 0.261
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Wind Burp 0.237
uva7b$keep_remove$keep
[1] "Defecate" "Itch"
$remove
[1] "Urinate" "Tickle"
uva7c <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva7cVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.457
Sneeze Cough 0.302
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.280
Urinate Defecate 0.254
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Hungry Thirsty 0.242
Vomit Sneeze 0.224
Heart Breathing 0.218
uva7c$keep_remove$keep
[1] "Cough" "Itch"
$remove
[1] "Sneeze" "Tickle"
uva8a <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva8aVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Urinate Defecate 0.237
Heart Breathing 0.218
Hungry Thirsty 0.213
uva8a$keep_removeNULL
uva8b <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva8bVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Urinate Defecate 0.277
Heart Breathing 0.267
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Wind Burp 0.206
uva8b$keep_removeNULL
uva9 <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva9Variable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.379
Wind Burp 0.373
Urinate Defecate 0.341
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Heart Breathing 0.285
Sneeze Cough 0.278
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Hungry Thirsty 0.24
uva9$keep_remove$keep
[1] "Defecate" "Wind" "Itch"
$remove
[1] "Urinate" "Burp" "Tickle"
uva10a <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva10aVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Tickle Itch 0.297
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Heart Breathing 0.209
uva10a$keep_removeNULL
uva10b <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva10bVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva10b$keep_removeNULL
uva10c <- EGAnet::UVA(data = data10c, cut.off = 0.3)Warning: Some variables did not belong to a dimension: Cough
Use caution: These variables have been removed from the TEFI calculation
uva10cVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva10c$keep_removeNULL
uva11a <- EGAnet::UVA(data = data11a, cut.off = 0.3)
uva11aVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Urinate Defecate 0.330
Heart Breathing 0.314
Wind Burp 0.302
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Muscles Pain 0.279
Tickle Itch 0.277
Hungry Thirsty 0.255
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Vomit Sneeze 0.201
uva11a$keep_remove$keep
[1] "Heart" "Defecate" "Burp"
$remove
[1] "Breathing" "Urinate" "Wind"
uva11b <- EGAnet::UVA(data = data11b, cut.off = 0.3)
uva11bVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Hungry Breathing 0.252
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Wind Burp 0.214
Temp Muscles 0.202
uva11b$keep_removeNULL
rez_uva <- c(uva1a$keep_remove$remove,
uva1b$keep_remove$remove,
uva2$keep_remove$remove,
uva3$keep_remove$remove,
uva4$keep_remove$remove,
uva5$keep_remove$remove,
uva6$keep_remove$remove,
uva7a$keep_remove$remove,
uva7b$keep_remove$remove,
uva7c$keep_remove$remove,
uva8a$keep_remove$remove,
uva8b$keep_remove$remove,
uva9$keep_remove$remove,
uva10a$keep_remove$remove,
uva10b$keep_remove$remove,
uva10c$keep_remove$remove,
uva11a$keep_remove$remove,
uva11b$keep_remove$remove)
sort(table(rez_uva), decreasing = TRUE) / length(rez_uva) rez_uva
Tickle Cough Urinate Breathing Burp Sneeze Wind
0.30434783 0.13043478 0.13043478 0.08695652 0.08695652 0.08695652 0.08695652
Defecate Itch
0.04347826 0.04347826
Discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant in UVA analysis.
data_all <- select(data_all, -Tickle)
data1a <- select(data1a, -Tickle)
data1b <- select(data1b, -Tickle)
data2 <- select(data2, -Tickle)
data3 <- select(data3, -Tickle)
data4 <- select(data4, -Tickle)
data5 <- select(data5, -Tickle)
data6 <- select(data6, -Tickle)
data7a <- select(data7a, -Tickle)
data7b <- select(data7b, -Tickle)
data7c <- select(data7c, -Tickle)
# data8a <- select(data8a, -Tickle)
# data8b <- select(data8b, -Tickle)
data9 <- select(data9, -Tickle)
data10a <- select(data10a, -Tickle)
data10b <- select(data10b, -Tickle)
data10c <- select(data10c, -Tickle)
data11a <- select(data11a, -Tickle)
data11b <- select(data11b, -Tickle)colors <- c("Heart" = "#F44336", "Breathing" = "#FF5722",
"Hungry" = "#FF9800", "Thirsty" = "#FFC107",
"Burp" = "#4CAF50", "Wind" = "#009688",
"Urinate" = "#63824C", "Defecate" = "#795548",
"Cough" = "#8BC34A", "Sneeze" = "#CDDC39",
"Bruise" = "#673AB7", "Blood Sugar" = "#3F51B5",
"Muscles" = "#2196F3", "Pain" = "#00BCD4",
"Sex arousal" = "#FF4081", "Temp" = "#9C27B0",
"Vomit" = "#76FF03", "Taste" = "#00E676",
"Affective touch" = "#FF1744", "Itch" = "#D500F9")make_hclust <- function(data) {
rez_pvclust <- pvclust::pvclust(data,
method.hclust = "complete",
method.dist = "correlation",
nboot = 1000, quiet=TRUE, parallel=TRUE)
dendrogram <- as.dist(1 - cor(data, use = "pairwise.complete.obs")) |>
hclust(method = "complete") |>
tidygraph::as_tbl_graph()
# Process Nodes
nodes <- as.list(dendrogram)$nodes |>
mutate(
Size = ifelse(label != "", 10, NA),
Item = str_replace(label, "_", " "),
idx = 1:nrow(as.list(dendrogram)$nodes))
# Central node
nodes[nodes$height == max(nodes$height), c("Item", "Size")] <- data.frame(Item="Central", Size=15)
# Process Edges
edges <- as.list(dendrogram)$edges
edges$linewidth = datawizard::rescale(nodes[edges$from, ]$height, to = c(0.1, 1))
p <- tidygraph::tbl_graph(nodes = nodes, edges = edges) |>
ggraph(layout = "dendrogram", circular = TRUE) +
# geom_edge_diagonal(strength = 0.7, linewidth = 1) +
geom_edge_elbow2(aes(edge_width=linewidth), color="#212121") +
geom_node_point(aes(filter=Item %in% c("Central"), size = Size), color="#212121") +
# geom_node_point(aes(filter=Group %in% c("Visceroception", "Awareness", "Deficit"), color=Group, size = Size)) +
geom_node_text(aes(
label = ifelse(Item != "", Item, NA),
x = x * 1.10,
y = y * 1.10,
filter = label != "",
angle = ifelse(
x >= 0,
asin(y) * 360 / 2 / pi,
360 - asin(y) * 360 / 2 / pi
),
hjust = ifelse(
x >= 0, 0, 1
))) +
geom_node_point(aes(filter = label != "", color=Item, size=Size), alpha = 1) +
# geom_node_text(aes(label=idx)) + # Debug
scale_edge_width_continuous(range=c(1, 3), guide = "none") +
scale_size_continuous(range=c(10, 10), guide = "none") +
scale_color_manual(values = colors,
breaks = names(colors)) +
ggtitle("Hierarchical Clustering Analysis (HCA)", subtitle = "Method = Correlation") +
coord_equal(clip = "off", xlim = c(-1.25, 1.25), ylim = c(-1.25, 1.25)) +
theme_void() +
# guides(color = guide_legend(override.aes = list(size = c(7.5, 5, 5, 3.5, 7.5, 5, 5, 3.5, 3.5, 3.5, 7.5, 5, 5, 3.5)))) +
theme(legend.text = ggtext::element_markdown(),
legend.title = element_blank(),
legend.position = "none",
# plot.title = element_blank(),
# plot.subtitle = element_blank()
plot.title = element_text(face="bold"))
list(pvclust = rez_pvclust, p = p)
}
rez_hclust <- make_hclust(data_all)
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)rez_hclust$pmake_ega <- function(data, cols = c("C1" = "#F44336", "C2" = "#FF9800", "C3" = "#795548", "C4" = "#CDDC39", "C5" = "#2196F3", "C6" = "#009688", "C7" = "#9C27B0")) {
rez <- EGAnet::bootEGA(
data = data,
seed = 123,
model = "glasso",
algorithm = "leiden",
EGA.type = "hierEGA",
type = "resampling",
plot.itemStability = FALSE,
verbose = FALSE,
allow.singleton = TRUE
)
table <- EGAnet::net.loads(rez$EGA$lower_order)$std |>
as.data.frame() |>
rownames_to_column("Item") |>
gt::gt() |>
gt::tab_header(title = "EGA Loadings") |>
gt::data_color(
columns = -Item,
method = "numeric",
colors = scales::col_numeric(
palette = c("red", "white", "green"),
domain = c(-1, 0, 1)
)) |>
gt::fmt_auto()
nodes <- rez$EGA$lower_order$dim.variables |>
rename(name = items) |>
mutate(dimension = paste0("C", dimension),
size = apply(EGAnet::net.loads(rez$EGA$lower_order)$std, 1, max)[name])
loadings <- rez$EGA$lower_order$network |>
as.data.frame() |>
rownames_to_column("to") |>
pivot_longer(-c(to), names_to = "from", values_to = "Loading") |>
filter(Loading > quantile(Loading, 1/3))
g <- tidygraph::tbl_graph(nodes = nodes, edges = loadings, directed = FALSE) |>
mutate(name = str_replace(name, "_", " "),
name = ifelse(name == "Sex arousal", "Sex", name))
set.seed(1)
layout <- ggraph::create_layout(g, layout = "fr", weights = abs(loadings$Loading))
xrange <- max(layout$x) - min(layout$x)
yrange <- max(layout$y) - min(layout$y)
xmin <- min(layout$x) - xrange * 0.05
xmax <- max(layout$x) + xrange * 0.05
ymin <- min(layout$y) - yrange * 0.05
ymax <- max(layout$y) + yrange * 0.05
p <- layout |>
ggraph() +
geom_edge_bend(aes(edge_width = Loading, edge_alpha = Loading), color = "#212121",
strength = 0.3) +
geom_node_point(aes(size = size, color = dimension), alpha = 0.95) +
geom_node_text(aes(label = name), size = 3, color = "white", fontface = "bold") +
scale_size_continuous(range = c(23, 28)) +
scale_edge_width(range = c(0.3, 4)) +
scale_edge_alpha(range = c(0.1, 0.9)) +
scale_color_manual(values = cols) +
labs(title = "Exploratory Graph Analysis (EGA)", subtitle = "Method = Leiden") +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(face="bold")) +
coord_cartesian(xlim =c(xmin, xmax), ylim = c(ymin, ymax))
list(rez = rez, table = table, p = p)
}
rez_ega <- make_ega(data_all, cols = c(
"C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]],
"C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]]))
stab <- EGAnet::itemStability(rez_ega$rez)rez_ega$table| EGA Loadings | |||||||
|---|---|---|---|---|---|---|---|
| Item | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Heart | 0.524 | 0.038 | 0.009 | 0.067 | 0.038 | 0.019 | 0.114 |
| Breathing | 0.524 | 0.193 | 0.046 | 0.185 | 0.055 | 0.002 | 0 |
| Hungry | 0.133 | 0.522 | 0.07 | 0.124 | 0 | 0.003 | 0.037 |
| Thirsty | 0.095 | 0.522 | 0.124 | 0.105 | 0.009 | 0.011 | 0.01 |
| Urinate | 0.021 | 0.179 | 0.563 | 0.15 | 0.038 | 0.015 | 0 |
| Defecate | 0.046 | 0.057 | 0.563 | 0.156 | 0.113 | 0.042 | 0 |
| Pain | 0.037 | 0.109 | 0.055 | 0.387 | 0.046 | 0 | 0.138 |
| Temp | 0.026 | 0.089 | 0.086 | 0.365 | 0.127 | 0.028 | 0.007 |
| Muscles | 0.103 | 0.024 | 0.039 | 0.359 | 0.058 | 0.089 | 0.127 |
| Sex_arousal | 0.073 | 0.065 | 0.071 | 0.298 | 0.079 | 0.078 | 0 |
| Affective_touch | 0.037 | 0 | 0.034 | 0.276 | 0.05 | 0 | 0.207 |
| Taste | 0.07 | 0.03 | 0.062 | 0.225 | 0.116 | 0.039 | 0.055 |
| Sneeze | 0.03 | 0 | 0.051 | 0.108 | 0.569 | 0.081 | 0.047 |
| Cough | 0.038 | 0.012 | 0.012 | 0.136 | 0.446 | 0.162 | 0.086 |
| Vomit | 0.047 | 0 | 0.094 | 0.188 | 0.293 | 0.043 | 0.042 |
| Wind | 0.008 | 0.009 | 0.055 | 0.116 | 0.104 | 0.572 | 0.059 |
| Burp | 0.019 | 0.01 | 0.005 | 0.099 | 0.184 | 0.572 | 0.111 |
| Bruise | 0.03 | 0 | 0 | 0.186 | 0.048 | 0.081 | 0.44 |
| Blood_Sugar | 0.067 | 0.049 | 0 | 0.077 | 0.02 | 0.025 | 0.387 |
| Itch | 0.019 | 0 | 0 | 0.138 | 0.076 | 0.033 | 0.298 |
rez_ega$pmake_efa <- function(data, n = 3, cols = c("F1" = "#009688", "F2" = "#795548", "F3" = "red")) {
n_fac <- n_factors(data)
rez_efa <- factor_analysis(data, n = n, sort = TRUE, rotation = "oblimin")
# Exctract loadings
loadings <- attributes(rez_efa)$loadings_long |>
select(from = Component, to=Variable, Loading) |>
mutate(Type = "Loading",
to = str_replace(to, "_", " "),
from = str_replace(from, "MR", "F"))
# Reorder factors
fct_order <- c("Pain", "Muscles", "Blood Sugar", "Bruise",
"Affective touch", "Sex arousal", "Temp", "Itch", "Tickle",
"Thirsty", "Hungry", "Breathing", "Heart", "Defecate", "Urinate",
"Sneeze", "Cough", "Wind", "Burp", "Vomit", "Taste")
loadings <- loadings |>
mutate(to = fct_relevel(to, fct_order[fct_order %in% unique(loadings$to)])) |>
arrange(to) |>
mutate(to = as.character(to))
# Get correlations between factors
phi <- attributes(rez_efa)$model$Phi
phi[upper.tri(phi, diag = TRUE)] <- NA
phi <- phi |>
as.data.frame() |>
rownames_to_column("from") %>%
pivot_longer(-from, names_to = "to", values_to = "Loading") |>
filter(!is.na(Loading)) |>
mutate(Type = "Correlation",
from = str_replace(from, "MR", "F"),
to = str_replace(to, "MR", "F"))
# Reverse so that arcs are on the same side
phi[paste(phi$from, phi$to) %in% c("F3 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F3")
phi[paste(phi$from, phi$to) %in% c("F2 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F2")
phi[paste(phi$from, phi$to) %in% c("F3 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F3")
phi[paste(phi$from, phi$to) %in% c("F4 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F4")
phi[paste(phi$from, phi$to) %in% c("F4 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F4")
phi[paste(phi$from, phi$to) %in% c("F4 F3"), c("from", "to")] <- data.frame(from = "F3", to = "F4")
loadings <- rbind(loadings, phi)
# Make layout
factor_names <- unique(loadings$from[loadings$Type == "Loading"])
item_names <- unique(loadings$to[loadings$Type == "Loading"])
nodes <- tibble(name = unique(c(loadings$from, loadings$to))) |>
mutate(type = case_when(
name %in% factor_names ~ "Factor",
name %in% item_names ~ "Item",
TRUE ~ "Other"
))
y_position <- seq(0.8, 0.2, length.out = length(factor_names))
layout <- nodes |>
mutate(x = ifelse(type == "Item", 0, 1),
y = case_when(
name == "F1" ~ y_position[1],
name == "F2" ~ y_position[2],
name == "F3" ~ y_position[3],
name == "F4" ~ y_position[4],
name == "F5" ~ y_position[5],
name == "F6" ~ y_position[6],
.default = NA
))
layout[layout$type == "Item", "y"] <- seq(0, 1, length.out = sum(layout$type == "Item"))
g <- tidygraph::as_tbl_graph(loadings, directed = TRUE) |>
left_join(layout, by = c("name"))
p_efa <- ggraph(g, layout = as.data.frame(layout)[c("x", "y")]) +
geom_edge_link(
aes(
edge_width = abs(Loading),
edge_alpha = abs(Loading),
# edge_color = Loading,
filter = Type == "Loading"
# label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
),
color = "#212121",
arrow = arrow(length = unit(2.5, "mm"), type = "closed"),
end_cap = circle(5, 'mm')
) +
geom_edge_arc(
aes(
edge_width = abs(Loading),
filter = Type == "Correlation",
label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
),
edge_color = "forestgreen",
angle_calc = "along",
label_dodge = unit(2.5, "mm"),
force_flip = TRUE,
strength = 0.2 # controls curvature
) +
geom_node_label(data=filter(layout, type == "Item"),
aes(fill = name),
label = paste(rep(" ", 20), collapse = ""),
size = 3,
hjust = 0.5,
nudge_x = -0.05,
label.r = unit(0.2, "lines"),
label.padding = unit(0.5, "lines"),
label.size = 0) +
geom_node_label(data=filter(layout, type == "Item"),
aes(label = name),
fill = NA,
size = 3,
hjust = 0.5,
vjust = 0.5,
color = "white",
fontface = "bold",
nudge_x = -0.05,
label.size = 0) +
geom_node_point(data=filter(layout, type == "Factor"),
aes(color = name),
size = 12) +
geom_node_text(data=filter(layout, type == "Factor"),
aes(label = name),
size = 4,
color = "white",
fontface = "bold") +
scale_edge_width(range = c(0.3, 2)) +
scale_edge_alpha(range = c(0.1, 1)) +
scale_fill_manual(values = colors) +
scale_color_manual(values = cols) +
labs(title = "Exploratory Factor Analysis (EFA)", subtitle = "Method = Oblimin") +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(face="bold")) +
coord_cartesian(xlim = c(-0.08, 1.08))
list(n_fac=n_fac, efa=rez_efa, p=p_efa)
}
rez_efa <- make_efa(data_all, n = 3, cols = c("F1" = colors[["Wind"]],
"F2" = colors[["Defecate"]],
"F3" = colors[["Heart"]]))Loading required namespace: GPArotation
plot(rez_efa$n_fac)rez_efa$pWe discarded:
data_allf <- select(data_all, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data1af <- select(data1a, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data1bf <- select(data1b, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data2f <- select(data2, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data3f <- select(data3, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data4f <- select(data4, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data5f <- select(data5, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data6f <- select(data6, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data7af <- select(data7a, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data7bf <- select(data7b, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data7cf <- select(data7c, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data8af <- select(data8a, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data8bf <- select(data8b, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data9f <- select(data9, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data10af <- select(data10a, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data10bf <- select(data10b, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data10cf <- select(data10c, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data11af <- select(data11a, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)
data11bf <- select(data11b, -Taste, -Affective_touch, -Vomit, -Itch,
-Sex_arousal, -Temp)rez_hclust <- make_hclust(data_allf)
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)rez_hclust$prez_ega <- make_ega(data_allf, cols = c(
"C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]],
"C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
stab <- EGAnet::itemStability(rez_ega$rez)rez_ega$table| EGA Loadings | |||||||
|---|---|---|---|---|---|---|---|
| Item | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Heart | 0.528 | 0.042 | 0.017 | 0.041 | 0.026 | 0.042 | 0.134 |
| Breathing | 0.528 | 0.203 | 0.065 | 0.054 | 0.01 | 0.17 | 0.01 |
| Hungry | 0.144 | 0.525 | 0.083 | 0.008 | 0.013 | 0.096 | 0.053 |
| Thirsty | 0.098 | 0.525 | 0.128 | 0.018 | 0.018 | 0.104 | 0.017 |
| Urinate | 0.032 | 0.192 | 0.572 | 0.043 | 0.025 | 0.116 | 0 |
| Defecate | 0.069 | 0.072 | 0.572 | 0.087 | 0.062 | 0.115 | 0 |
| Sneeze | 0.057 | 0.009 | 0.089 | 0.563 | 0.104 | 0.087 | 0.03 |
| Cough | 0.054 | 0.021 | 0.036 | 0.563 | 0.176 | 0.12 | 0.078 |
| Wind | 0.014 | 0.024 | 0.071 | 0.117 | 0.576 | 0.074 | 0.067 |
| Burp | 0.032 | 0.015 | 0.016 | 0.179 | 0.576 | 0.098 | 0.127 |
| Muscles | 0.125 | 0.035 | 0.067 | 0.082 | 0.107 | 0.494 | 0.16 |
| Pain | 0.054 | 0.135 | 0.09 | 0.066 | 0.008 | 0.494 | 0.198 |
| Bruise | 0.042 | 0 | 0 | 0.051 | 0.092 | 0.281 | 0.488 |
| Blood_Sugar | 0.075 | 0.057 | 0 | 0.023 | 0.035 | 0.065 | 0.488 |
rez_ega$prez_efa <- make_efa(data_allf, n = 3, cols =
c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]],
"F3" = colors[["Heart"]]))
plot(rez_efa$n_fac)rez_efa$pmake_cfa <- function(data) {
# Main models
m_g <- '
IAS =~ Burp + Wind + Cough + Sneeze + Hungry + Thirsty + Urinate + Defecate + Heart + Breathing + Pain + Muscles
'
m_ega <- '
CoughSneeze =~ Cough + Sneeze
UrinateDefecate =~ Urinate + Defecate
BruiseBloodsugar =~ Bruise + Blood_Sugar
WindBurp =~ Wind + Burp
HeartBreath =~ Heart + Breathing
MusclesPain =~ Muscles + Pain
HungryThirsty =~ Hungry + Thirsty
'
m_efa <- '
UrinateDefecate =~ Urinate + Defecate
Sudden =~ Cough + Sneeze + Wind + Burp
Homeostasis =~ Heart + Breathing + Muscles + Pain + Hungry + Thirsty + Bruise + Blood_Sugar
'
m_hclust <- '
HeartBreath =~ Heart + Breathing
BruiseBloodsugar =~ Bruise + Blood_Sugar
HungryThirsty =~ Hungry + Thirsty
Sudden =~ Cough + Sneeze + Wind + Burp
Unpleasant =~ Urinate + Defecate + Muscles + Pain
'
# Test bifactor models
m_egag <- paste0(m_ega, "\nIAS =~ CoughSneeze + WindBurp + BruiseBloodsugar + UrinateDefecate + HeartBreath + MusclesPain")
m_efag <- paste0(m_efa, "\nIAS =~ UrinateDefecate + Sudden + Homeostasis")
m_hclustg <- paste0(m_hclust, "\nIAS =~ HeartBreath + BruiseBloodsugar + HungryThirsty + Sudden + Unpleasant")
fit_g <- lavaan::cfa(m_g, data = data, std.lv = TRUE)
fit_ega <- lavaan::cfa(m_ega, data = data, std.lv = TRUE)
fit_egag <- lavaan::cfa(m_egag, data = data, std.lv = TRUE)
fit_efa <- lavaan::cfa(m_efa, data = data, std.lv = TRUE)
fit_efag <- lavaan::cfa(m_efag, data = data, std.lv = TRUE)
fit_hclust <- lavaan::cfa(m_hclust, data = data, std.lv = TRUE)
fit_hclustg <- lavaan::cfa(m_hclustg, data = data, std.lv = TRUE)
# Test
bf <- test_bf(fit_ega, fit_g, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg) |>
as.data.frame() |>
select(Name = Model, BF)
# test <- lavaan::lavTestLRT(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg)
test <- compare_performance(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg,
metrics = c("BIC", "RMSEA", "Chi2", "CFI")) |>
select(-BIC_wt, -Model) |>
merge(bf, by = "Name") |>
arrange(BIC) |>
gt::gt() |>
gt::tab_header(title = "CFA Model Comparison") |>
gt::data_color(
columns = c(BIC, Chi2, RMSEA),
method = "numeric",
palette = c("green", "white")) |>
gt::data_color(
columns = c(CFI),
method = "numeric",
palette = c("white", "green")) |>
gt::data_color(
columns = c(BF),
method = "numeric",
palette = c("red", "yellow", "white", "yellow", "green"),
domain = c(0, 0.3, 1, 3, 10000000)) |>
gt::fmt_number(columns = c(RMSEA, Chi2, CFI), decimals = 3)
param <- parameters::parameters(fit_ega)
perf <- performance::performance(fit_ega)
# Graph
nodes <- tidySEM::get_nodes(fit_ega) |>
mutate(name = str_replace(name, "_", " "))
edges <- tidySEM::get_edges(fit_ega) |>
mutate(from = str_replace(from, "_", " "),
to = str_replace(to, "_", " "),
est_std = as.numeric(est_std)) |>
filter(from != to)
g <- tidygraph::tbl_graph(nodes = nodes, edges = edges, directed = TRUE)
set.seed(1)
layout <- ggraph::create_layout(g, layout = "fr", weights = est_std)
p_cfa <- ggraph(layout) +
geom_edge_bend(aes(edge_width = abs(est_std), edge_alpha = abs(est_std),
label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
filter = arrow == "none"),
color = "forestgreen",
angle_calc = "along",
check_overlap = TRUE,
label_dodge = unit(0, "mm"),
label_size = unit(3, "mm"),
strength = 0.3) +
geom_edge_link(aes(edge_width = abs(est_std), edge_color = abs(est_std),
label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
filter = arrow != "none"),
angle_calc = "along",
label_dodge = unit(3, "mm"),
arrow = arrow(length = unit(2.5, "mm"), type = "open", ends = "first"),
lineend = "round",
linejoin = "bevel",
start_cap = circle(9, 'mm')) +
geom_node_point(aes(color = name, shape=shape), alpha = 0.97, size = 20) +
geom_node_text(data = filter(layout, shape == "rect"), aes(label = name),
size = 3, color = "white", fontface = "bold") +
scale_size_continuous(range = c(10, 20)) +
scale_edge_width(range = c(0.3, 4)) +
scale_edge_alpha(range = c(0.1, 0.9)) +
scale_edge_color_gradient(low = "#EEEEEE", high = "#212121") +
scale_color_manual(values = c(colors,
"CoughSneeze" = "#CDDC39",
"UrinateDefecate" = "#795548", "BruiseBloodsugar" = "#673AB7",
"WindBurp" = "#009688", "HeartBreath" = "#F44336",
"MusclesPain" = "#2196F3", "HungryThirsty" = "#FF9800")) +
scale_shape_manual(values = c("oval" = 16, "rect" = 15)) +
labs(title = "Confirmatory Factor Analysis (CFA)", subtitle = "Method = Maximum Likelihood") +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(face="bold"))
p_cfa
list(test = test, param = param, perf = perf, p=p_cfa)
}
rez_cfa0 <- make_cfa(data_allf)Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: When comparing models, please note that probably not all models were fit
from same data.
When comparing models, please note that probably not all models were fit
from same data.
rez_cfa0$test| CFA Model Comparison | |||||
|---|---|---|---|---|---|
| Name | BIC | RMSEA | Chi2 | CFI | BF |
| fit_ega | -158722.7 | 0.036 | 2,290.274 | 0.984 | NA |
| fit_egag | -154419.8 | 0.055 | 6,738.245 | 0.952 | 0 |
| fit_hclust | -148000.0 | 0.079 | 13,126.977 | 0.905 | 0 |
| fit_hclustg | -146829.7 | 0.079 | 14,349.064 | 0.896 | 0 |
| fit_g | -145655.5 | 0.125 | 26,541.561 | 0.783 | 0 |
| fit_efa | -144933.2 | 0.083 | 16,266.333 | 0.883 | 0 |
| fit_efag | -144933.2 | 0.083 | 16,266.333 | 0.883 | 0 |
rez_cfa0$perf# Indices of model performance
Chi2(56) | p (Chi2) | Baseline(91) | p (Baseline) | GFI | AGFI | NFI
-------------------------------------------------------------------------
2290.274 | < .001 | 1.380e+05 | < .001 | 0.990 | 0.981 | 0.983
Chi2(56) | NNFI | CFI | RMSEA | RMSEA CI | p (RMSEA) | RMR | SRMR
-----------------------------------------------------------------------------
2290.274 | 0.974 | 0.984 | 0.036 | [0.034, 0.037] | > .999 | 0.001 | 0.020
Chi2(56) | RFI | PNFI | IFI | RNI | Loglikelihood | AIC
---------------------------------------------------------------------
2290.274 | 0.973 | 0.605 | 0.984 | 0.984 | 79615.220 | -1.591e+05
Chi2(56) | BIC | BIC_adjusted
------------------------------------
2290.274 | -1.587e+05 | -1.589e+05
rez_cfa0$param# Loading
Link | Coefficient | SE | 95% CI
-----------------------------------------------------------------------
CoughSneeze =~ Cough | 0.17 | 1.16e-03 | [0.17, 0.17]
CoughSneeze =~ Sneeze | 0.15 | 1.10e-03 | [0.15, 0.15]
UrinateDefecate =~ Urinate | 0.17 | 1.23e-03 | [0.17, 0.17]
UrinateDefecate =~ Defecate | 0.16 | 1.19e-03 | [0.16, 0.16]
BruiseBloodsugar =~ Bruise | 0.22 | 2.11e-03 | [0.21, 0.22]
BruiseBloodsugar =~ Blood_Sugar | 0.16 | 1.94e-03 | [0.15, 0.16]
WindBurp =~ Wind | 0.18 | 1.29e-03 | [0.18, 0.18]
WindBurp =~ Burp | 0.20 | 1.32e-03 | [0.20, 0.20]
HeartBreath =~ Heart | 0.14 | 1.49e-03 | [0.14, 0.14]
HeartBreath =~ Breathing | 0.17 | 1.44e-03 | [0.17, 0.18]
MusclesPain =~ Muscles | 0.15 | 1.19e-03 | [0.15, 0.15]
MusclesPain =~ Pain | 0.15 | 1.23e-03 | [0.14, 0.15]
HungryThirsty =~ Hungry | 0.18 | 1.62e-03 | [0.18, 0.18]
HungryThirsty =~ Thirsty | 0.18 | 1.58e-03 | [0.18, 0.19]
Link | z | p
-------------------------------------------------
CoughSneeze =~ Cough | 147.35 | < .001
CoughSneeze =~ Sneeze | 136.55 | < .001
UrinateDefecate =~ Urinate | 135.80 | < .001
UrinateDefecate =~ Defecate | 136.25 | < .001
BruiseBloodsugar =~ Bruise | 102.29 | < .001
BruiseBloodsugar =~ Blood_Sugar | 80.90 | < .001
WindBurp =~ Wind | 141.15 | < .001
WindBurp =~ Burp | 150.94 | < .001
HeartBreath =~ Heart | 95.60 | < .001
HeartBreath =~ Breathing | 119.51 | < .001
MusclesPain =~ Muscles | 124.42 | < .001
MusclesPain =~ Pain | 119.27 | < .001
HungryThirsty =~ Hungry | 112.05 | < .001
HungryThirsty =~ Thirsty | 115.23 | < .001
# Correlation
Link | Coefficient | SE | 95% CI
---------------------------------------------------------------------------
CoughSneeze ~~ UrinateDefecate | 0.56 | 5.94e-03 | [0.55, 0.57]
CoughSneeze ~~ BruiseBloodsugar | 0.55 | 7.32e-03 | [0.53, 0.56]
CoughSneeze ~~ WindBurp | 0.72 | 4.80e-03 | [0.71, 0.73]
CoughSneeze ~~ HeartBreath | 0.52 | 6.73e-03 | [0.51, 0.53]
CoughSneeze ~~ MusclesPain | 0.67 | 5.99e-03 | [0.65, 0.68]
CoughSneeze ~~ HungryThirsty | 0.46 | 7.02e-03 | [0.45, 0.47]
UrinateDefecate ~~ BruiseBloodsugar | 0.40 | 7.83e-03 | [0.38, 0.41]
UrinateDefecate ~~ WindBurp | 0.51 | 6.10e-03 | [0.50, 0.52]
UrinateDefecate ~~ HeartBreath | 0.54 | 6.71e-03 | [0.52, 0.55]
UrinateDefecate ~~ MusclesPain | 0.65 | 6.14e-03 | [0.64, 0.66]
UrinateDefecate ~~ HungryThirsty | 0.65 | 6.15e-03 | [0.64, 0.66]
BruiseBloodsugar ~~ WindBurp | 0.58 | 7.09e-03 | [0.57, 0.59]
BruiseBloodsugar ~~ HeartBreath | 0.50 | 8.08e-03 | [0.48, 0.51]
BruiseBloodsugar ~~ MusclesPain | 0.74 | 7.37e-03 | [0.73, 0.76]
BruiseBloodsugar ~~ HungryThirsty | 0.45 | 8.30e-03 | [0.43, 0.47]
WindBurp ~~ HeartBreath | 0.46 | 6.86e-03 | [0.44, 0.47]
WindBurp ~~ MusclesPain | 0.63 | 6.03e-03 | [0.62, 0.65]
WindBurp ~~ HungryThirsty | 0.43 | 7.01e-03 | [0.42, 0.45]
HeartBreath ~~ MusclesPain | 0.64 | 6.88e-03 | [0.63, 0.66]
HeartBreath ~~ HungryThirsty | 0.64 | 6.90e-03 | [0.63, 0.66]
MusclesPain ~~ HungryThirsty | 0.63 | 6.90e-03 | [0.62, 0.65]
Link | z | p
-----------------------------------------------------
CoughSneeze ~~ UrinateDefecate | 94.23 | < .001
CoughSneeze ~~ BruiseBloodsugar | 74.98 | < .001
CoughSneeze ~~ WindBurp | 150.87 | < .001
CoughSneeze ~~ HeartBreath | 77.12 | < .001
CoughSneeze ~~ MusclesPain | 111.01 | < .001
CoughSneeze ~~ HungryThirsty | 65.50 | < .001
UrinateDefecate ~~ BruiseBloodsugar | 51.02 | < .001
UrinateDefecate ~~ WindBurp | 83.50 | < .001
UrinateDefecate ~~ HeartBreath | 79.79 | < .001
UrinateDefecate ~~ MusclesPain | 105.72 | < .001
UrinateDefecate ~~ HungryThirsty | 105.96 | < .001
BruiseBloodsugar ~~ WindBurp | 81.97 | < .001
BruiseBloodsugar ~~ HeartBreath | 61.56 | < .001
BruiseBloodsugar ~~ MusclesPain | 100.55 | < .001
BruiseBloodsugar ~~ HungryThirsty | 54.31 | < .001
WindBurp ~~ HeartBreath | 66.62 | < .001
WindBurp ~~ MusclesPain | 105.09 | < .001
WindBurp ~~ HungryThirsty | 61.95 | < .001
HeartBreath ~~ MusclesPain | 93.45 | < .001
HeartBreath ~~ HungryThirsty | 93.47 | < .001
MusclesPain ~~ HungryThirsty | 91.62 | < .001
rez_cfa0$pThe structural analysis seem to converge on the idea of small clusters (of pairs or triplets) that are potentially inter-related parts of larger clusters (although with unstable associations). Best way to assess the values of this new granular structure is to assess its predictive value against convergent measures, and see if it’s superior to a unique score (-> Study 2).
fig2 <- (rez_hclust$p | rez_ega$p) / (rez_efa$p | rez_cfa0$p)
ggsave("figures/Figure2.png", fig2, width=14, height=14, dpi=300, bg="white")rez_hclust <- make_hclust(data1a)
rez_ega <- make_ega(data1a, cols = c(
"C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]],
"C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1a, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]],
"F3" = colors[["Heart"]], "F4" = colors[["Itch"]]))
rez_hclust$p / rez_ega$p / rez_efa$pplot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)stab <- EGAnet::itemStability(rez_ega$rez)rez_ega$table| EGA Loadings | ||||||||
|---|---|---|---|---|---|---|---|---|
| Item | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| Heart | 0.504 | 0.125 | 0.046 | 0 | 0.093 | 0.114 | 0 | 0 |
| Breathing | 0.504 | 0.14 | 0.032 | 0.098 | 0.158 | 0.123 | 0 | 0.036 |
| Hungry | 0.2 | 0.511 | 0.12 | 0.004 | 0.028 | 0.064 | 0 | 0.034 |
| Thirsty | 0.074 | 0.511 | 0.166 | 0 | 0.068 | 0.13 | 0.005 | 0.053 |
| Defecate | 0.044 | 0.093 | 0.632 | 0.107 | 0.077 | 0.036 | 0.059 | −0.045 |
| Urinate | 0.026 | 0.253 | 0.421 | 0.035 | 0.167 | 0.035 | 0.046 | −0.02 |
| Taste | 0.032 | 0.018 | 0.21 | 0.105 | 0.18 | 0.149 | 0.082 | 0.028 |
| Sneeze | 0 | 0.006 | 0.078 | 0.676 | 0.094 | 0.001 | 0.092 | 0.044 |
| Cough | 0.053 | 0 | 0.014 | 0.377 | 0.11 | 0 | 0.171 | 0.121 |
| Vomit | 0.088 | 0 | 0.177 | 0.299 | 0.066 | 0.067 | 0 | 0.054 |
| Muscles | 0.084 | 0.022 | 0.13 | 0.03 | 0.454 | 0.032 | 0.029 | 0.093 |
| Pain | 0.126 | 0.032 | 0.133 | 0.067 | 0.438 | 0.195 | 0 | 0.143 |
| Temp | 0.086 | 0.056 | 0.117 | 0.125 | 0.264 | 0.245 | 0 | 0.038 |
| Sex_arousal | 0.103 | 0.042 | 0.055 | 0.006 | 0.069 | 0.435 | 0.119 | 0.007 |
| Affective_touch | 0.059 | 0.086 | 0.058 | 0.027 | 0.205 | 0.435 | 0.011 | 0.152 |
| Wind | 0 | 0.007 | 0.189 | 0.051 | 0 | 0.261 | 0.571 | 0.058 |
| Burp | 0 | 0 | 0.007 | 0.203 | 0.034 | 0 | 0.571 | 0.105 |
| Blood_Sugar | 0.013 | 0 | −0.02 | 0.04 | 0.033 | 0.033 | 0.015 | 0.455 |
| Bruise | 0 | 0.108 | −0.044 | 0.08 | 0.193 | 0.044 | 0.068 | 0.401 |
| Itch | 0.034 | 0 | 0.027 | 0.076 | 0.074 | 0.224 | 0.07 | 0.388 |
plot(rez_efa$n_fac)rez_hclust <- make_hclust(data1af)
rez_ega <- make_ega(data1af, cols = c(
"C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]],
"C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1af, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]],
"F3" = colors[["Heart"]], "F4" = colors[["Itch"]]))
rez_cfa1a <- make_cfa(data1af)
rez_hclust$p / rez_ega$p / rez_efa$pplot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)stab <- EGAnet::itemStability(rez_ega$rez)rez_ega$table| EGA Loadings | |||||||
|---|---|---|---|---|---|---|---|
| Item | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Heart | 0.511 | 0.128 | 0.036 | 0.009 | 0.01 | 0.108 | 0 |
| Breathing | 0.511 | 0.16 | 0.048 | 0.068 | 0.015 | 0.151 | 0.045 |
| Hungry | 0.202 | 0.514 | 0.105 | 0.015 | 0 | 0.042 | 0.048 |
| Thirsty | 0.091 | 0.514 | 0.178 | 0 | 0.018 | 0.041 | 0.089 |
| Urinate | 0.029 | 0.253 | 0.565 | 0.038 | 0.05 | 0.185 | −0.017 |
| Defecate | 0.08 | 0.108 | 0.565 | 0.098 | 0.087 | 0.09 | −0.042 |
| Sneeze | 0.033 | 0.02 | 0.12 | 0.574 | 0.117 | 0.091 | 0 |
| Cough | 0.071 | 0 | 0.021 | 0.574 | 0.171 | 0.06 | 0.157 |
| Wind | 0.021 | 0.024 | 0.136 | 0.071 | 0.577 | 0.028 | 0.089 |
| Burp | 0.012 | 0 | 0.008 | 0.22 | 0.577 | 0.041 | 0.089 |
| Muscles | 0.091 | 0.042 | 0.118 | 0.058 | 0.04 | 0.523 | 0.125 |
| Pain | 0.183 | 0.044 | 0.107 | 0.061 | 0.014 | 0.523 | 0.18 |
| Bruise | 0.01 | 0.128 | −0.031 | 0.059 | 0.093 | 0.222 | 0.502 |
| Blood_Sugar | 0.033 | 0 | −0.013 | 0.052 | 0.032 | 0.053 | 0.502 |
plot(rez_efa$n_fac)rez_cfa1a$test| CFA Model Comparison | |||||
|---|---|---|---|---|---|
| Name | BIC | RMSEA | Chi2 | CFI | BF |
| fit_ega | -1635.181 | 0.026 | 72.825 | 0.988 | NA |
| fit_egag | -1613.162 | 0.059 | 180.403 | 0.924 | 1.654961e-05 |
| fit_hclustg | -1594.961 | 0.065 | 210.828 | 0.905 | 1.846821e-09 |
| fit_hclust | -1592.739 | 0.062 | 182.492 | 0.921 | 6.080789e-10 |
| fit_efag | -1564.540 | 0.073 | 253.472 | 0.877 | 4.576430e-16 |
| fit_efa | -1564.540 | 0.073 | 253.472 | 0.877 | 4.576430e-16 |
| fit_g | -1538.631 | 0.106 | 325.704 | 0.784 | 1.082552e-21 |
TODO.
rbind(
mutate(rez_cfa0$perf, Sample = "All Samples"),
mutate(rez_cfa1a$perf, Sample = "Sample 1a")
# TODO
) |>
select(Sample, Chi2, RMSEA, CFI, "SRMR") # Indices of model performance
Sample | Chi2 | RMSEA | CFI | SRMR
----------------------------------------------
All Samples | 2290.274 | 0.036 | 0.984 | 0.020
Sample 1a | 72.825 | 0.026 | 0.988 | 0.029
Note the usage of descriptive factor names relating directly to the items to avoid abstraction.
# Add empirical variables
add_empirical <- function(data, sample = "Sample 1a") {
x <- data.frame(
Original = rowMeans(data),
HungryThirsty = (data$Hungry + data$Thirsty) / 2,
MusclesPain = (data$Muscles + data$Pain) / 2,
WindBurp = (data$Wind + data$Burp) / 2,
UrinateDefecate = (data$Urinate + data$Defecate) / 2,
BreathingHeart = (data$Heart + data$Breathing) / 2)
if("Blood_Sugar" %in% names(data)) {
x$BruiseBlood <- (data$Blood_Sugar + data$Bruise) / 2
} else {
x$BruiseBlood <- data$Bruise
}
if("Cough" %in% names(data)) {
x$CoughSneeze <- (data$Cough + data$Sneeze) / 2
} else {
x$CoughSneeze <- data$Sneeze
}
x$Sample <- sample
x
}
scores <- list(
sample1a = add_empirical(data1a, "Sample 1a"),
sample1b = add_empirical(data1b, "Sample 1b"),
sample2 = add_empirical(data2, "Sample 2"),
sample3 = add_empirical(data3, "Sample 3"),
sample4 = add_empirical(data4, "Sample 4"),
sample5 = add_empirical(data5, "Sample 5"),
sample6 = add_empirical(data6, "Sample 6"),
sample7a = add_empirical(data7a, "Sample 7a"),
sample7b = add_empirical(data7b, "Sample 7b"),
sample7c = add_empirical(data7c, "Sample 7c"),
sample8a = add_empirical(data8a, "Sample 8a"),
sample8b = add_empirical(data8b, "Sample 8b"),
sample9 = add_empirical(data9, "Sample 9"),
sample10a = add_empirical(data10a, "Sample 10a"),
sample10b = add_empirical(data10b, "Sample 10b"),
sample10c = add_empirical(data10c, "Sample 10c"),
sample11a = add_empirical(data11a, "Sample 11a"),
sample11b = add_empirical(data11b, "Sample 11b")
)
save(scores, file = "../data/scores.RData")